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Summary 


A  conventional  frequency  domain  helicopter-borne  electro¬ 
magnetic  (HEM)  system  can  be  used  to  map  sea  ice  keels  with  a  rea¬ 
sonable  degree  of  accuracy.  A  preliminary  interpretation  of  the  ac¬ 
quired  data  can  be  made  manually  with  the  help  of  a  nomogram  or 
automated  with  the  use  of  a  table  look-up  routine  on  a  small  com¬ 
puter.  Such  data  may  also  be  more  accurately  interpreted  with  the 
use  of  an  adaptation  of  Occam's  inversion.  This  scheme  allows  for  the 
practical  non-uniqueness  of  the  inverse  solution  but  selects  the 
smoothest  keel  shape  that  is  consistent  with  the  field  data.  The 
inversion  method  is  much  more  computationally  intensive  than  the 
table  look-up  technique.  While  the  latter  can  be  implemented  on  a 
small  computer  to  form  an  interactive  in-flight  interpretation 
system,  the  inversion  technique  involves  many  forward 
computations  and,  for  the  present,  is  best  reserved  for  past  flight 
data  analysis.  It  is  possible  that  this  difficulty  can  be  resolved  with 
the  use  of  specialized  computing  equipment. 

In  the  strict  sense  both  the  proposed  interpretation  techniques 
are  only  suitable  for  use  on  data  acquired  over  two  dimensional 
features  whose  strike  length  (measured  in  a  direction  perpendicular 
to  the  flight  line)  is  much  greater  than  the  flight  height.  An  examina¬ 
tion  of  the  anomalies  for  three-dimensional  keels  however,  reveals 
that  good  data  interpretation  is  possible  whenever  the  keel  strike 
length  exceeds  the  system  height  by  a  factor  of  three. 
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Airborne  Electromagnetic  Mapping  of  Sea  Ice  Keels 

A  lex  Becker  and  Guimin  Liu 


1.0  Introduction 

Sea  ice  is  a  unique  environment  encountered  in  most  Arctic 
work.  This  includes  the  transportation  of  vehicles  through  ice- 
covered  waters,  the  construction  of  offshore  drilling  structures,  and 
the  safe  operation  of  submarines.  In  such  circumstances,  knowledge 
of  the  thickness  and  properties  of  sea  ice  is  important.  This  is  true 
especially  for  ice  keels,  which  constitute  a  local  downward 
indentation  of  the  ice-water  interface. 

Recently,  Canpolar  consultants  (1985)  reviewed  the  possible 
techniques  for  remotely  measuring  sea  ice  thickness,  among  these 
they  included  airborne  impulse  radar  and  airborne  electromagnetics. 
The  impulse  radar  technique  has  been  very  successful  in  some  areas 
(Kovacs  and  Morey,  1985).  It  may  however,  occasionally  yield 
erroneous  data  when  saline  moisture  zones  exist  within  the  ice 
cover.  Becker  et  al.  (1983),  on  the  other  hand,  examined  the 
feasibility  of  using  frequency-domain  airborne  electromagnetics  to 
determine  sea  ice  thickness.  At  low  frequencies,  sea  ice  is 
practically  transparent  to  electromagnetic  waves  and  the  observed 
secondary  magnetic  field  can  be  used  to  estimate  the  distance  from 
the  EM  system  boom  to  the  ice-water  interface.  At  the  same  time,  a 
laser  altimeter  installed  on  the  boom  measures  the  distance  to  the 
surface  of  the  ice,  so  that  the  difference  of  these  two  distances  is 
the  sea  ice  thickness.  This  technique  was  successfully  used  in 
mapping  the  average  sea  ice  thickness  in  Prudhoe  Bay,  Alaska,  but  it 
failed  in  an  area  of  a  multi-year  pressure  keel  because  of  the 
inappropriate  one-dimensional  (1-D)  techniques  used  to  interpret 
data  (Kovacs,  et  al.,  1986;  Becker,  et  al.,  1987).  In  order  to  recover 


.  e  geometry  of  the  keel,  two-  or  three-dimensional  interpretation 
techniques  are  required. 

In  this  report  we  primarily  concern  ourselves  with  the 
development  of  techniques  for  interpreting  AEM  data  for  two- 
dimensional  keels  which  are  assumed  to  be  infinitely  long  in  a 
direction  perpendicular  to  the  flight  line.  The  computational 
methods  used  to  construct  the  necessary  forward  solutions  are 
outlined  in  the  Appendix.  These  serve  either  for  the  construction  of 
interpretation  nomograms  or  in  the  imputation  of  a  numerical 
inversion  scheme.  In  order  to  establish  the  validity  of  the  two- 
dimensional  interpretation  scheme  we  also  examine  the  effects  of 
finite  keel  length  on  the  observable  anomalies. 

In  terms  of  computer  time,  the  modeling  technique  we  employ 
is  particularly  well  suited  for  the  interpretation  of  the  AEM  data 
collected  over  sea  ice  keels.  Using  the  conventional  finite  element 
method,  30  minutes  CRAY  II  CPU  time  is  required  to  compute  the  AEM 
system  response  along  a  traverse  line  over  an  ice  keel.  However,  for 
a  similar  problem,  using  an  approach  rooted  in  potential  theory,  the 
computation  takes  about  10  seconds  on  IBM  3090  (equivalent  to 
about  3  seconds  on  CRAY)  with  a  gain  more  than  two  orders  of 
magnitude  in  speed.  We  will  show  that  the  assumption  of  a  perfectly 
conducting  sea  water  model  is  not  a  significant  d.awback  since  it 
can  be  used  for  any  system  operation  frequency  greater  than  30  KHz. 

2.0  AEM  Anomalies  for  Ice  Keels 

Let  us  first  briefly  examine  the  airborne  electromagnetic 
response  to  two-dimensional  sea  ice  keels  as  a  function  of  their 
size  and  shape  (cf.  Fig.1).  For  these  calculations,  the 
electromagnetic  system  has  a  coil  separation  of  6.5m  and  is  "flown" 
25m  above  the  upper  ice  surface  which  is  flat.  With  the  exception  of 
the  zone  containing  the  keel,  the  sea  ice  is  5m  thick  and  is  assumed 
to  have  negligible  electrical  conductivity.  The  data  is  calculated  for 
both  the  vertical-coaxial  coil  configuration  (HX  system)  and  the 
horizontal-coplanar  one  (HZ  system).  The  keel  is  assumed  to  have 
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tne  shape  of  a  Gaussian  curve  and  its  parameters  and  geometry  are 
shown  in  Figure  2. 

Figure  3  shows  the  HX  system  response  for  two  different 
keels  with  parameters  A  =  12m,  and  W  =  28m  (solid  line)  or  14m 
(dashed  line).  The  observed  secondary  magnetic  field  is  expressed  in 
parts  per  million  (ppm)  of  the  primary  field  at  the  receiver.  There  is 
a  big  anomaly  in  the  system  response  due  to  either  of  these  ice 
keels.  For  a  keel  width  W  =  28m,  the  maximum  anomaly  is  about  35% 
of  the  background  response  and  we  define  this  to  be  the  "percent 
anomaly".  In  order  to  characterize  the  shape  of  the  observed 
anomaly,  we  define  the  "anomaly  width",  q,  to  be  the  width  of  the 
anomaly  at  one  half  of  its  maximum  value.  For  a  keel  width  W  =  28m, 
the  anomaly  width  q  is  45m.  When  the  keel  width  W  is  halved  and 
other  conditions  are  kept  unchanged,  the  HX  system  anoma'y  is 
significantly  reduced.  In  this  case,  the  maximum  anomaly  (or  percent 
anomaly)  decreases  by  37%,  and  the  anomaly  width  decreases  by 
33%. 

The  situation  for  the  HZ  system  is  similar.  Figure  4  shows  the 
HZ  system  response  due  to  the  same  two  ice  keels  that  we-e 
previously  examined.  For  a  keel  defined  by  A  =  12m  and  W  =28m,  the 
percent  anomaly  for  the  HZ  system  response  is  32%.  The  anomaly 
width  q  is  66m  in  this  case,  which  is  much  larger  than  that  seen  for 
the  HX  system.  When  the  keel  width  is  halved,  the  HZ  system 
anomaly  also  decreases  greatly.  In  this  case,  the  anomaly  amplitude 
(or  percent  anomaly)  drops  by  32%  and  the  anomaly  width  decreases 
by  20%.  From  the  comparison  of  Figures  3  and  4,  we  can  see  that  the 
HX  system  response  tracks  the  shape  of  the  the  keel  closely,  while 
the  HZ  system  anomaly  is  much  wider. 

Next,  the  keel  drawdown  A  is  halved  while  the  keel  width  is 
kept  constant  at  W  =  28m.  The  corresponding  HX  system  response  is 
shown  in  Figure  5  (dashed  line),  and  HZ  system  response  is  shown  in 
Figure  6  (dashed  line).  The  anomaly  amplitude  contracts  by  29%  for 
the  HX  system  and  by  37%  for  the  HZ  system.  However,  the 
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corresponding  decrease  in  the  anomaly  width  is  small,  8%  for  the  HX 
system,  and  3%  for  the  HZ  system. 

From  the  above  model  studies,  we  find  that  the  amplitude  of 
the  anomaly  (or  accordingly  the  percent  anomaly)  is  sensitive  to 
both  the  thickness  and  the  width  of  the  keel.  Thus,  it  is  essentially  a 
function  of  the  area  of  the  cross  section  of  the  keel.  In  contrast,  the 
anomaly  width  is  primarily  related  to  the  keel  width.  It  is  not 
sensitive  to  the  keel  thickness  as  long  as  the  shape  of  the  keel  does 
not  change.  It  is  also  noticed  that  the  HX  system  response  is  more 
sensitive  to  keel  shape  than  the  HZ  system  response.  The  curve  of  HX 
system  response  resembles  the  keel  shape  and  the  change  in  the 
anomaly  width  due  to  the  change  of  the  keel  width  is  much  larger 
than  that  for  the  HZ  system  (33%  compared  to  20%). 

It  is  worthwhile  to  compare  the  result  calculated  with  the 
above  technique  to  that  obtained  by  the  finite  element  method  (Lee 
and  Morrison,  1985)  that  allows  a  finite  conductivity  for  both  the 
sea  ice  and  the  sea  water  as  well  as  a  low  operation  frequency.  This 
is  done  in  Figure  7  for  a  triangular  sea  ice  keel.  The  ice  is  uniformly 
5m  thick  except  for  the  keel  area  and  it  has  a  conductivity  of  0.002 
S/m.  The  conductivity  of  the  sea  water  is  4  S/m  and  the  HX  system 
flies  at  30m  above  the  surface  of  the  ice  and  operates  at  2500  Hz. 
The  dashed  line  in  Figure  7  (upper  part)  represents  the  in-phase 
component  of  the  response,  and  the  solid  line  corresponds  to  the 
system  response  that  would  be  observed  if  the  sea  water 
conductivity  were  infinite  and  sea  ice  conductivity  were  zero. 
Notice  the  similarity  of  these  two  curves  and  that  the  percent 
anomalies  and  anomaly  widths  are  almost  identical  for  these  two 
curves.  Thus  it  appears  appropriate  to  use  the  perfect  conductor 
model  to  interpret  low  frequency  data  observed  in  practice. 

3.0  Charts  for  Data  Interpretation 

From  the  above  analysis  of  numerical  data,  it  >s  seen  that  the 
percent  anomaly  is  a  function  of  the  cross  area  of  a  sea  ice  keel  and 
the  anomaly  w;dth  is  primarily  a  function  of  the  keel  width.  To 
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interpret  field  data,  we  need  a  strategy  to  relate  the  observed 
electromagnetic  anomaly  to  the  keel  geometry.  In  terms  of  the 
smooth  keel  used  in  this  analysis  (see  Figure  2),  we  need  to 
estimate  the  two  keel  parameters  A  and  W  from  the  observed 
anomaly  of  the  system  response.  In  the  following,  we  are  going  to 
concentrate  on  the  HX  system  configuration  which  appears  superior 
in  terms  of  sensitivity  of  the  keel  geometry.  Parallel  analysis  can 
be  carried  out  for  the  HZ  system.  . 

From  the  computation  and  analysis  of  the  HX  system  response 
for  model  keels  with  systematic  values  of  A  and  W,  one  can 
construct  an  interpretation  chart  as  shown  in  Figure  8.  The  vertical 
axis  of  the  chart  is  the  percent  anomaly,  while  the  horizontal  axis 
represents  the  ratio  of  the  anomaly  width  q  to  the  flight  height  h 
above  the  flat  part  of  the  ice-water  interface.  Here,  the  two  sets  of 
parametric  curves  intersect  each  other  clearly  and  are  well 
separated.  The  solid  lines  are  for  constant  values  of  the  drawdown 
A,  and  the  dashed  lines  are  for  constant  keel  width  W.  The  charted 
values  of  a  and  w  are  the  keel  drawdown  "A"  and  width  "W" 
normalized  by  h,  the  system  height  above  the  flat  seawater  surface. 
As  shown,  the  percent  anomaly  decreases  with  the  decrease  of  the 
keel  drawdown  and  keel  width.  In  fact,  the  line  a  =  0.8  tends  to  zero 
fast  and  intersects  the  lines  a  =  0.4,  a  *  0.2,  etc..  For  purpose  of 
clarity,  this  is  not  shown  in  Figure  8.  Nonetheless,  this  does  not 
pose  a  serious  problem  in  practice  since  such  narrow  and  sharp  keels 
are  highly  improbable.  Although  the  interpretation  chart  is 
constructed  for  h  =  30m,  it  can  be  used  for  the  range  h  =  25m  -  50m 
with  an  error  less  than  ±  4%.  . 

As  pointed  out  at  the  end  of  last  section,  the  anomaly  shape 
for  low  frequency  data  is  almost  identical  to  that  obtained  in  case 
of  the  perfect  conductor  model.  This  makes  the  interpretation  chart 
shown  in  Figure  8  applicable  for  a  wide  range  of  frequencies.  We 
expect  that  the  interpretation  chart  is  useful  for  the  frequency 
range  between  1-100  kHz.  A  similar  interpretation  chart  can  be 
constructed  for  the  HZ  system  response.  It  is  shown  in  Figure  9. 


Now  consider  the  2500  Hz  theoretical  data  shown  in  Figure  7 
(dashed  line)  for  an  assymmetric  triangular  keel.  From  the  in-phase 
component  of  the  AEM  response,  the  percent  anomaly  and  anomaly 
width  are  calculated  to  be  7%  and  27m  respectively.  Thus  q/h  = 
27/35  =  0.72.  We  draw  this  point  in  the  interpretation  chart  (point  C 
in  Figure  8)  and  find  the  corresponding  values  of  a  and  w  to  be  0.1 
and  0.35  respectively.  Hence  A  =  a  x  h  =  3.5m,  and  W  =  wxh»  12.3m. 
The  estimated  keel  is  drawn  (symmetrically  about  the  maximum 
anomaly  of  the  data)  in  the  lower  portion  of  the  illustration  (dashed 
line  in  Figure  7).  We  can  see  that  the  size  of  the  interpreted  keel  is 
close  to  that  of  the  model  keel.  But  since  we  have  assumed  that  the 
keel  is  symmetric  in  the  interpretation,  the  position  of  the 
maximum  of  the  interpreted  keel  is  misplaced  3.5m  to  the  right. 

The  field  data  collected  over  an  ice  keel  in  the  Prudhoe  Bay  by 
Geotech  Ltd.  for  CRREL  (Kovacs  et  al.,  1986)  are  interpreted  next. 
The  AEM  system  used  consists  of  two  pairs  of  vertical  coaxial  coils 
(HX  system)  and  two  pairs  of  horizontal  coplanar  coils  (HZ  system). 
The  former  operate  at  930  Hz  and  4158  Hz  respectively,  and  the 
latter  operate  at  530  Hz  and  16290  Hz  respectively.  The 
transmitter-receiver  separation  of  each  coil  pair  is  about  6.5m. 

A  part  of  the  930  Hz  and  16290  Hz  data  for  line  F6L3  is  shown 
in  Figure  10.  The  altitude  is  the  distance  from  the  system  boom  to 
the  ice  surface  measured  by  a  laser  altimeter.  Note  that  the 
quadrature  component  of  the  16290  Hz  data  is  of  bad  quality.  As  we 
can  see  in  the  figure,  the  data  are  highly  correlated  with  the 
altitude  as  expected.  First,  we  interpreted  the  data  using  a  1-D 
technique  (Becker  et  al,  1987).  The  result  shows  an  average  of  3m 
thick  ice  (see  Figure  11)  but  no  ice  keel  is  apparent  in  the 
interpretation.  However,  we  notice  that  there  is  an  anomaly  in  the 
system  response  from  fiducial  numbers  2668  -  2675,  which  can  not 
be  related  to  the  small  altitude  variation  in  that  area.  Moreover,  the 
1-D  interpreted  results  also  show  thicker  ice  in  that  region.  For  the 
930  Hz  data,  the  anomaly  width  q  is  found  to  be  31.4m  and  the 
percent  anomaly  is  6.5%.  Since  h  is  about  38.5m  (altitude  +  average 
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ice  thickness),  q/h  =  0.83.  The  corresponding  point  is  found  in  the 
interpretation  chart  (point  D  in  Figure  8),  which  gives  a  «=  0.08  and  w 
=  0.42.  Hence  A  =  a  x  h  =  3.08m,  W  =  w  x  h  =  16.2m.  The  interpreted 
keel  is  plotted  symmetrically  about  the  maximum  anomaly  in  Figure 
11  (dashed  line).  The  solid  line  in  the  figure  is  the  average  of  the 
drill  hole  measurements  along  three  parallel  lines  11.5  meters 
apart.  As  we  can  see,  the  interpreted  keel  is  an  good  approximation 
of  the  real  feature.  Note,  however,  that  the  effect  of  the  variation  of 
the  altitude  of  the  system  boom  has  been  neglected  in  the 
interpretation. 

The  interpretation  of  4158  Hz  data  yields  similar  result  to  the 
above.  But  difficulty  was  encountered  in  interpreting  the  HZ  system 
data  which  were  acquired  at  530  and  16290  Hz.  As  shown  in  Figure 
10,  there  is  also  an  anomaly  in  the  16290  Hz  data  from  fiducial 
numbers  2668  -  2675.  If  we  calculate  the  percent  anomaly  and 
anomaly  width  in  this  case,  the  corresponding  point  will  fall  off  the 
interpretation  chart  shown  in  Figure  9  and  can  not  be  interpreted. 

Generally  speaking,  the  chart  interpretation  method  gives  an 
overall  estimate  of  the  size  and  shape  of  an  ice  keel.  In  order  to 
recover  the  detail  shape  of  an  ice  keel,  we  consider  in  the  next 
section  an  iterative  technique  for  data  inversion.  The  chart 
interpretation  can  be  used  to  obtain  an  initial  model  for  the 
iterative  technique  where  the  variation  of  the  system  altitude  can 
be  easily  accounted  for. 

4.0  inversion  of  AEM  Data 

Now  consider  the  determination  of  sea  ice  keel  shape  by  data 
inversion.  In  order  to  overcome  the  problem  associated  with  1-D 
data  interpretation,  we  present  a  two-dimensional  (2-D)  AEM  data 
inversion  scheme  where  ice  thickness  may  vary  along  the  flight 
direction  of  the  AEM  system.  In  the  perpendicular  direction, 
however,  the  geometry  of  the  ice/water  interface  is  assumed  to  be 
invariant.  Although  the  actual  ice  bottom  topography  is  three 
dimensional  (3-D),  most  ice  keels  do  exhibit  a  dominant  extension  in 


8 


one  direction  which  is  called  the  strike  because  they  are  formed  by 
the  interaction  of  two  ice  plates  (Canpolar,  1985).  Indeed,  as  will  be 
shown  in  the  next  section,  if  the  strike  length  of  an  ice  keel  of 
Gaussian  shape  is  greater  than  three  times  the  AEM  system  flight 
height  2-D  interpretation  may  be  quite  appropriate. 

The  interpretation  scheme  that  we  have  developed  is  based  on 
the  principle  of  the  Occam's  inversion  presented  by  Constable  et  al. 
(1987).  Due  to  practical  non-uniqueness  of  the  inverse  problem, 
there  may  be  many  sea  ice  keels  that  may  fit  the  observed  data 

within  a  specified  error.  This  scheme  yields  the  smoothest  among 
these.  The  smooth  keel  shows  the  basic  features  of  the  actual  keel 

although  small  bumps  on  its  surface  may  not  appear  in  the  inversion 
results.  Because  the  process  of  the  data  acquisition  however 

constitutes  a  low-pass  filter,  it  is  also  possible  that  the  information 
needed  to  define  the  keel  in  detail  is  absent  from  the  acquired  data. 

We  first  examine  the  theory  of  constrained  smooth  inversion 
and  then  use  it  to  interpret  synthetic  field  data  which  are  generated 
with  our  numerical  modeling  program  (cf.  Appendix).  Both  noise  free 
and  artificially  contaminated  data  are  used  to  test  the  scheme. 

Finally  we  propose  a  procedure  to  apply  this  inversion  scheme  to  low 
frequency  data.  In  addition  to  this  field  data  sets  (collected  in 
Prudhoe  Bay  by  Geotech, Ltd.  in  1985)  are  also  interpreted  using  this 
procedure. 


4.1  Theory 

Let  x  denote  the  traverse  direction  of  the  AEM  system,  y  the 
strike  direction  of  the  ice  keel  and  z  the  downward  pointing  vertical. 
The  interface  of  the  ice  and  the  water  is  described  by  z  =  h(x).  Our 
purpose  here  is  to  reconstruct  this  interface  from  the  airborne 
electromagnetic  data  d(x)  while  the  upper  surface  cf  the  ice  is 
mapped  by  a  laser  device.  The  difference  in  elevation  between  the 
bottom  and  the  top  surface  of  the  ice  is  the  required  ice  thickness. 
The  data  d(x)  can  be  either  the  horizontal  or/and  the  vertical 
secondary  magnetic  field  recorded  during  the  flight. 
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Now  consider  h(x)  as  the  only  unknown.  This  is  true  when  the 
two  assumptions  used  in  our  previous  modeling  algorithm  (see  the 
Appendix)  are  valid.  These  are:  1)  the  sea  ice  is  transparent  to  the 
electromagnetic  wave,  and  2)  the  sea  water  can  be  treated  as  a 
perfect  electric  conductor.  With  these  considerations,  the  AEM  data 
can  be  written  as 

d(x)  =  F[h(x)]  ...  (1) 

where  F[  ]  is  the  functional  for  forward  modeling.  Equation  (1) 
immediately  reveals  that  the  inversion  problem  is  actually  one¬ 
dimensional  (1-D)  with  the  direction  of  variation  along  x  instead  of 

along  z.  Fortunately  there  is  a  large  amount  of  geophysical  literature 
dealing  with  1-D  inversion.  One  particular  scheme  that  well  suits  our 
needs  is  Occam’s  inversion  presented  by  Constable  et  al.  (1987). 

The  basis  for  Occam's  inversion  is  the  search  for  the  smoothest 
solution  that  fits  the  observed  data  within  a  specified  tolerance.  This 
principle  applies  particularly  well  to  our  work  because 
electromagnetic  wave  propagation  in  sea  water  is  a  diffusion 

process,  so  that  the  resolution  of  sharp  edges  on  the  ice/water 
interface  can  not  be  expected  from  the  data.  Furthermore,  any 
information  on  the  sharp  edges  is  contained  in  the  high  wave- 

number  range  of  the  data,  which  is  contaminated  by  noise.  Because 
inversion  implies  downward  continuation  (Parker,  1977),  the 
attempts  to  reconstruct  the  fine  structures  of  the  interface  will 
amplify  any  noise  and  result  in  unstable  solutions.  Hence  a  stable 
solution  is  necessarily  smooth. 

In  order  to  find  the  smoothest  solution,  let  us  first  define  the 
roughness  of  the  ice/water  interface.  Physically  the  rougher  the 

interface  is,  the  larger  the  magnitude  of  its  derivatives.  Thus  we 
define 
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to  be  the  roughness  of  the  interface  h(x);  (pi ,  p2)  define  the  lateral 
keel  extent  outside  of  which  the  interface  is  assumed  to  be  flat.  In 
the  discrete  sense, 

N 

R  =  Z(hrh1-i)2  =  H3h||2 

-  ...  (2) 


where  hi  =  h(xi  )  and  II  II  denotes  the  /  2  norm  ,  and 
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Note  that  pi=  xi,  p2  =  xn  and  that  the  points  xi  ,  X2  ,  ....  xn  are  usually 
equally  spaced. 


Suppose  that  there  are  M  data  points  recorded  over  an  ice  keel, 
dj  =  d(Xj),  j=l,  2,  ...  M.  The  corresponding  computed  predictions  from 
the  discrete  model  h  are  Fj(h).  The  goodness  of  fit  of  the  predictions 
to  the  actual  data  can  be  evaluated  using  the  least-squares  criterion 


H  =  £[dj-FJ(h)]2=||d-F(h)||: 


i-i 


...  (3) 


where 


d  =  (d i,  dj  ...,  df,i)T 

F(h)  =  (F1(h),F2(h),...,FM(h)  / 

The  ice/water  interface  h  can  only  vary  within  some  physical 
bounds.  If  we  assume  that  the  z  =  o  plane  is  chosen  such  that  it 
coincides  with  the  flat  part  of  that  interface  then  a  reasonable  lower 
bound  is: 


h,  >  0, 


i  = 


...  (4) 


since  the  ice  keel  protrudes  downward.  An  upper  bound  can  also  be 
set  for  each  individual  case 


hj  <  Tj ,  i  =  1, 2, ....  N  ...  (5) 

In  the  Arctic,  small  ice  keels  may  protrude  several  meters  into  the 
water,  whereas  large  keels  can  protrude  tens  of  meters  (Lowry  and 
Wadhams,  1979).  The  values  of  Ti  should  be  estimated  for  each 
specific  ice  keel  encountered.  Note  that  h]  =  hN  =0  should  be 
included  in  the  constraints  in  order  for  the  solution  to  be  smooth  at 
the  two  end  points.  This  condition  can  be  met  by  simply  letting  Tj  = 
T\  =  0. 

Now  the  mathematical  problem  to  be  solved  can  be  stated  as 
follows:  find  a  solution  h  which  minimizes  the  roughness  R  and 
brings  the  misfit  E  within  an  acceptable  tolerance,  while  the  bound 
constraints  (4)  and  (5)  are  satisfied.  Without  the  bound  conditions  (4) 
and  (5),  this  problem  is  exactly  identical  to  the  one  solved  by 
Constable  et  al.  (1987). 

The  condition  for  the  data  misfit  is 

l|d  .  F(h)]|2 <  E,  (6) 

where  E.  is  the  tolerance.  If  we  treat  this  inequality  as  an  equality 
and  apply  the  method  of  Lagrange  multipliers,  the  above  problem 
can  be  reduced  to  the  minimization  of 

U  =  ||3h||2+  p-1  ( ||d  -  F(h)!}2  —  E.)  _  (7) 

with  constraints  (4)  and  (5).  Here  p_1  is  the  Lagrange  multiplier.  As 
interpreted  by  Constable  et  al,  p  is  a  smoothing  parameter.  The  larger 
the  (i  is,  the  less  the  solution  is  affected  by  the  misfit.  On  the 
contrary,  if  p  is  small,  the  data  misfit  is  minimized  with  little 
influence  from  the  roughness  term. 

The  original  problem  can  now  be  solved  with  the  following 
procedures:  Solve  the  above  minimization  problem  for  a  series  of  p 


values  to  obtain  a  set  of  solutions  of  the  ice/water  interface.  Among 
these  solutions,  choose  the  one  which  satisfies  the  tolerance  condition 
(6).  If  more  than  one  solution  satisfies  (6),  choose  the  one  with  the 
largest  p  value  for  this  corresponds  to  the  smoothest  solution 
desired. 

Such  solutions  can  not  however  be  easily  obtained  by  the  direct 
minimization  of  the  objective  function  U  (  equation  (7)  )  since  the 
minimization  is  non-linear.  It  is  first  necessary  to  transform  the  non¬ 
linear  problem  into  a  problem  of  quadratic  programming,  for  which 
existing  mathematical  tools  can  be  used. 

Let  us  linearize  the  F(h)  about  an  initial  model  hO 

F(h)  =  F(h°)  +  J  A  ...  (g) 

Here  h  =  h°  +  A,  and  the  Jacobian  matrix  is 

3F !  BFi 

3h  j  3h 2  dh^ 

5F2  3F2  3F2 

3h  j  dh~2  '  ’  5h  >sf 

3FM3FM  3Fm 
3h  j  3hY  ’  ’  ’  3h^ 

Substituting  (8)  into  (7)  and  dropping  the  constant  term  F  E*  ,  we 
obtain 

U  =  ||3h|j2-i-  p_1  |  jd  .  J  h  1 1 2 

where  H  =  d  -  F(h°)  +  J  h°  js  the  modified  data.  Rearrangement  of 
the  above  equation  gives 

V  =  IpU-dTd  =  IhTH  h  -  CTh 


Here  V  is  the  new  objective  function  and 
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h  =  paTa  +  jtj 

C  =  JTd 

Note  that  H  is  a  symmetric  positive -definite  matri  x.  The 

minimization  of  the  new  objective  function  V  is  equivalent  to  the 
minimization  of  U  for  a  fixed  value  of  p  . 

Now  that  the  new  objective  function  is  in  quadratic  form  the 

problem  of  optimization  with  bound  constraints  (4)  and  (5)  can  be 
solved  using  quadratic  programming  (Gill  et  al.,  1981).  There  are 
subroutines  available  in  existing  mathematical  software  libraries 
which  can  be  used  for  this  purpose.  These  are  E04AF  in  the  NAG 
Fortran  Library  and  VE04A  in  the  Harwell  Fortran  Library.  We 
selected  VE04A  because  of  its  simplicity. 

The  smoothest  solution  can  now  be  actually  obtained  in  the 
following  way.  Starting  from  an  initial  model  hO,  solve  the 

minimization  problem  for  different  p  values.  From  these  solutions 
choose  the  one  that  minimizes  the  actual  misfit  E  instead  of  V. 
(Minimizing  V  may  result  in  divergence  of  the  solution  (Constable  et 
al.,  1987)).  Use  this  solution  for  ihe  next  initial  model  and  iterate 
until  the  solution  for  the  ice/water  interface  that  brings  the  misfit 

below  a  specified  tolerance  is  found. 

The  initial  model  can  be  chosen  arbitrarily  since  it  does  not 
affect  the  convergence  of  the  inversion  to  the  smoothest  solution. 
This  is  one  of  the  beauties  of  the  smooth  inversion  scheme  which  sets 
out  to  seek  a  unique  solution.  In  our  problem  we  set  the  initial  model 
to  be  a  flat  ice  cake,  i.e.,  =  (0,  0 . 0)T  . 


4.2  Inversion  of  synthetic  data 

In  this  section  we  will  apply  the  above  scheme  to  invert  some 
synthetic  data  which  are  generated  with  the  fast  modeling  algorithm 
(see  Appendix).  To  test  the  stability  of  the  inversion  scheme  white 


noise  will  also  be  added  to  the  numerical  data.  Most  inversion 
schemes  require  the  knowledge  of  solutions  to  the  forward  problem 
in  the  computation  of  the  Jacobian  matrix.  Here  these  are  obtained 
with  the  fast  forward  modeling  algorithm.  The  partial  derivatives  in 
the  Jacobian  are  computed  using  the  forward  finite-difference  of  two 
numerical  solutions.  Thus  one  iteration  of  the  inversion  process 
requires  N+l  forward  calculations  for  M  data  points. 

Before  considering  the  inversion  results  let  us  first  describe  the 
geometry  of  the  problem  for  all  the  synthetic  models.  Unless 
otherwise  specified,  the  transmitter  and  the  receiver  are  both  small 
circular  loops  separated  by  6.5  meters,  and  both  are  maintained  at 
25  meters  above  the  upper  ice  surface.  For  convenience  the  vertical 
co-axial  coil  system  is  referred  to  as  the  HX  system  because  the  axes 
of  the  coils  are  in  the  x  direction  while  the  horizontal  co-planar  coil 
system  is  referred  to  as  the  HZ  system.  The  receiver  measures  the 
secondary  magnetic  field  which  is  expressed  in  parts  per  million  of 

the  received  primary  field.  The  inversion  is  performed  for  the  HX 
and  the  HZ  system  data  independently  although  the  joint  inversion 
can  be  easily  accomplished.  Except  in  the  keel  area  the  ice  thickness 
is  assumed  to  be  five  meters.  Note  that  for  the  synthetic  data 

generation  we  assume  the  induction  limit  to  hold  so  that  the  system 
frequency  is  not  involved.  Applications  to  low  frequency  data, 
however,  will  be  shown  in  the  next  section. 

The  position  of  the  end  points  pi  and  P2  of  the  ice  keel  must 
now  be  chosen.  These  can  be  arbitrary  as  long  as  they  are  located 
outside  the  range  of  the  keel.  They  should  not  be  too  far  apart 
however,  because  the  computation  of  the  Jacobian  matrix  can  be 
quite  expensive  and  convergence  of  the  inversion  may  be  slow.  From 

our  experience,  locating  these  two  points  at  the  one  quarter  of  the 

peak  HX  system  anomaly  points,  and  at  the  two  points  at  half  of  the 
peak  HZ  system  anomaly  yields  good  results  .  Note  that  this  condition 
may  need  to  be  relaxed  for  field  data  because  they  may  be  acquired 
at  fluctuating  flight  heights,  which  will  distort  the  shape  of  the 
observed  anomaly. 
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The  sampling  interval  of  h(x),  Ax,  is  taken  to  he  3  meters  for  all 
the  examples  shown  in  this  report.  The  data  d(x)  is  sampled  at  an 

interval  of  3.5  meters  unless  otherwise  specified,  which  at  a 
helicopter  speed  of  68  knots  corresponds  to  a  recording  rate  of  10 

samples/second. 

Figure  12(a)  shows  the  inversion  results  of  the  HX  system  data 
for  a  triangular  model  keel  (Model  1),  which  is  symmetric  and  is 

located  between  60  and  90  meters  along  the  profile  The  keel 
drawdown  at  the  peak  is  5  meters  and  its  two  sides  slope  at  IS 

degrees  from  the  horizontal.  The  solid  line  in  the  upper  graph  is  the 
synthetic  data  corresponding  to  the  ice  keel  shown  in  solid  line  in  the 
lower  graph,  while  the  dashed  line  is  the  system  response  from  the 
inverted  keel  shown  in  dashed  line  in  the  lower  graph.  Similar 
results  for  the  HZ  system  data  are  shown  in  Figure  12(b).  In  the 
inversion  the  two  end  points  pi  and  P2  are  taken  at  45  and  105 
meterc  r-’rr° rtively .  The  tolerance  criterion  for  convergence  is  set  at 
1  ppm  for  the  HX  system  data  and  2  ppm  for  the  HZ  system  data. 
Three  iterations  in  the  inversion  yield  the  convergent  solutions 
shown  in  Figure  12  from  an  initial  guess  of  a  flat  ice  cake  that  is  5 
meters  thick.  As  we  can  see  the  interpreted  keels  have  smooth 

vertices  as  can  be  expected  for  the  constrained  smooth  inversion.  The 
interpreted  keels  are  about  1.5  meters  too  shallow,  and  the 
independent  results  for  the  HX  and  the  HZ  system  data  are  almost 
identical. 

In  the  next  example  we  invert  the  synthetic  data  for  an 

irregular  trapezoidal  ice  keel  (Model  2).  The  steep  side  of  the  keel  is 
34  degrees  from  the  horizontal  while  the  slope  of  the  other  side  is  16 
degrees  (Figure  13).  The  keel  protrudes  6  meters  dowm  into  the 
water  with  a  flat  bottom  of  12  meters.  Again  both  the  HX  and  the  HZ 
system  data  are  inverted  independently  and  the  results  shown  were 
obtained  after  three  iterations.  Misfit  for  the  HX  system  data  is  again 

held  below  lppm,  and  that  for  the  HZ  system  data  is  below  2ppm. 

For  the  HX  system  the  inverted  maximum  keel  thickness  is  correct 
although  the  two  vertices  on  the  steep  side  are  smoothed  out.  The 


other  side  however  is  imaged  correctly.  For  the  HZ  system,  the 
results  are  similar  except  thm  the  maximum  keel  thickness  is  a  half 
meter  too  small. 

Now  white  noise  of  5%  to  10%  of  the  peak  anomaly  is  added  to 
the  Model  1  and  Model  2  theoretical  data.  Here  system  noise  may  be 
small  but  "geological  noise”  is  estimated  to  be  about  5  percent  of  the 
anomaly  amplitude.  The  tolerance  criterion  for  the  inversion  is  set  at 
the  RMS  noise  level  and  the  convergence  is  usually  achieved  within  4 
-5  iterations.  The  inversion  results  are  shown  in  Figures  14  -  16.  As 
demonstrated  in  these  figures  the  inversion  still  yields  good  results 
except  at  the  10%  noise  level  where  the  keel  interpreted  from  the  HZ 
system  data  is  flat  at  the  bottom  (Figure  15).  Our  experience  with 
noisy  data  shows  that  the  inversion  results  obtained  from  the  HX 
system  data  are  usually  better  than  those  obtained  from  the  HZ 
system  data.  In  practice  5%  noise  is  probably  excessive. 

We  next  examine  the  case  where  the  height  of  the  system 
varies  during  a  flight,  as  is  usually  the  case  in  practice.  The  model 
keel  is  a  symmetric  trapezoid  in  shape  with  a  maximum  drawdown 
of  3  meters  (Model  3).  The  flat  bottom  is  12  meters  wide  and  the  two 
sides  of  the  keel  slope  at  12  degrees  from  the  horizontal  (Figure  17). 
The  system  height  over  the  flat  part  of  the  ice/water  interface, 
shown  in  Figure  17(a),  varies  between  37  and  41  meters.  This  height 
variation  was  taken  from  a  test  flight  in  Prudhoe  Bay.  Here  it  is 
assumed  that  the  transmitter  and  the  receiver  are  always  at  the 
same  level  although  in  practice  the  instrument  pod  may  tilt  in  space. 
The  inversion  results  for  the  HX  and  the  HZ  system  response  are 
shown  in  Figures  17(b)  and  (c).  As  can  be  seen  the  inverted  keels 
agree  very  well  with  the  model  keel  in  this  case.  The  final  data  misfit 
is  0.3ppm  for  the  HX  system  and  0.8ppm  for  the  HZ  system. 


4.3  Application  to  low  frequency  data 


The  inversion  procedure  can  also  be  extended  to  the 
interpretation  of  low  frequency  field  data.  As  mentioned  earlier,  the 
synthetic  data  used  above  were  generated  in  the  induction  limit 
where  the  secondary  field  only  exhibits  an  in-phase  component.  In 
practice  however,  this  can  not  be  realized  because  the  sea  water  is 
not  a  superconductor  and  the  operating  frequency  of  the  system 
must  be  low  enough  for  the  EM  wave  to  penetrate  the  ice  freely. 
Hence  the  measured  secondary  magnetic  field  has  both  an  in-phase 
and  a  quadrature  component.  Ideally  the  in-phase  and  quadrature 
components  can  be  inverted  simultaneously  using  the  above  scheme 
if  a  fast  modeling  algorithm  for  a  general  conductivity  distribution 
exists.  The  finite-element  and  finite-difference  methods  have  been 
very  successful  in  solving  such  problems  (Lee  and  Morrison,  1 9 S 5 : 
Stover  and  Greenfield,  1983).  The  computational  costs  associated 
with  these  methods  however,  are  prohibitive  so  that  they  can  not  be 
used  in  solving  this  problem  as  too  many  forward  computations  must 
to  be  earned  out  even  for  one  iteration. 

Although  the  fast  modeling  algorithm  that  we  developed 
generates  theoretical  data  at  the  induction  limit  it  can  be  used  to 

invert  low-frequency  data.  For  data  collected  at  30  kHz  the 
arithmetic  sum  of  the  in-phase  and  quadrature  components  gives  a 
good  approximation  of  the  data  that  would  be  obtained  at  the 

induction  limit  (Becker  et  al,  1983).  This  sum  can  be  directly 

interpreted  using  the  previous  scheme,  but  it  must  first  be 
multiplied  by  a  scale  factor  prior  to  interpretation.  The  same  scheme 
can  be  used  since  the  anomaly  shape  changes  but  little  with 

frequency  (Becker  and  Liu,  1987).  The  scale  factor  can  be  computed 
as  the  ratio  of  the  response  at  the  induction  limit  to  the  sum  of  the 
in-phase  and  quadrature  components  at  that  frequency  for  a  1-D 
model.  Note  that  this  factor  varies  with  the  system  altitude.  If  the 
altitude  does  not  change  much  in  a  flight  over  the  ice  keel  this  factor 
can  be  taken  as  a  constant. 


Let  us  now  consider  an  example  of  synthetic  data  of  the  HX 
system  at  the  frequency  of  2500  Hz  which  was  generated  by  the 
finite-element  method  (Lee  and  Momson,  1985).  The  conductivities 
of  the  ice  and  water  are  0.002  S/m  and  4  S/m  respectively.  The  ice  is 
uniformly  5  meters  thick  except  at  the  keel  area  where  it  protrudes 
5  meters  into  the  sea  water.  The  shape  of  the  keel  is  step-triangular 
as  shown  in  Figure  18  and  the  HX  system  is  flown  at  30  meters 
above  the  ice  surface.  The  scaled  sum  of  the  in-phase  and  quadrature 
components  of  the  theoretical  da'a  and  the  inversion  results  are  also 
shown  in  Figure  18.  There  are  ten  unevenly  spaced  data  points  and 
the  scale  factor  is  1.044  in  this  case.  As  we  can  see  the  inverted  ice 
keel  is  close  to  the  model  although  it  is  somewhat  shallower.  We 
consider  this  an  encouraging  success  of  the  application  of  the  above 
proposed  procedure. 

We  have  not  yet  had  the  opportunity  to  experiment  with  large 
quantities  of  low  frequency  data.  Thus  it  is  not  clear  in  which 
frequency  range  the  above  procedure  can  be  used  to  yield 
reasonable  results.  Furthermore  the  role  of  the  water  conductivity 
has  not  yet  been  investigated  although  we  suspect  tt'.at  its  main 
effect  is  to  decrease  the  resolution  of  the  inversion  results  expected 
from  the  superconductor  assumption.  It  is  worthwhile  however  to 
attempt  an  inversion  of  the  data  collected  in  the  Prudhoe  Bay  by 
Geotech  Ltd.  in  1985  (Kovacs,  et  al.,  1987).  The  data  is  for  line  F6L3 
collected  with  the  HZ  system  at  the  frequency  of  16290  Hz.  Because 
the  quadrature  component  is  of  very  poor  quality  only  the  in-phase 
component  was  scaled  by  a  factor  of  1.104  and  interpreted.  The 
system  altitude,  the  scaled  in-phase  component  and  the  inversion 
results  are  all  shown  in  Figure  19.  The  inversion  was  performed  with 
the  two  end  points  of  the  ice  keel  fixed  at  pi  =  54  meters  and  p; 
=  102  meters.  The  background  ice  thickness  was  fixed  at  3.3  meters 
which  is  the  average  from  the  1-D  interpretation. 

The  interpreted  ice  thickness  is  1.5  meters  too  shallow  as 
compared  to  the  drill-hole  measurements  (solid  line  in  the  bottom  of 
Figure  19)  in  the  keel  area.  However  the  keel  is  clearly  visible  in  the 


inversion  result,  which  is  encouraging  since  it  is  much  better  than 
the  1-D  inversion  that  shows  but  little  of  the  keel  (Becker  ci  al., 
1987).  The  inversion  was  stopped  at  an  RMS  data  misfit  of  15  ppm 
since  more  iterations  could  not  reduce  the  misfit.  There  may  be 
several  causes  for  such  a  large  misfit: 

1)  The  effects  of  the  pitch  and  roll  of  the  system  during  the  flight 
was  not  accounted  for  in  this  inversion. 

2)  Th3  time  constant  of  the  instrument  may  be  too  large  for  such  a 
small  keel  to  be  fully  represented  in  the  data.  The  inversion  treats 
the  data  as  being  recorded  with  a  zero  time  constant,  which  ignores 
the  integral  effect  due  to  the  instrument.  This  effect  has  been 
carefully  studied  by  Becker  and  Cheng  (1987). 

3)  The  top  surface  of  the  ice  was  treated  as  a  flat  plane  and  the 
effect  of  the  associated  ridge  was  not  considered. 

4)  The  assumed  2-D  structure  of  the  ice  may  be  incorrect. 

5)  For  this  frequency  channel  the  operational  system  noise  w-as 
large.  The  quadrature  component  was  quite  erratic  and  was  not  used 
in  the  inversion. 

Errors  due  to  1)  can  be  diminished  by  improving  the  computer 
code.  Errors  due  to  2)  and  5)  can  be  reduced  to  a  large  extent  by 

improving  the  instrument.  But  errors  due  to  4)  can  not  be  reduced  in 
the  present  2-D  inversion  scheme.  To  do  this  a  3-D  inversion 

algorithm  needs  to  be  developed  which  can  be  an  extension  of  the 
present  scheme.  The  process  will  be  computer  intensive  and 

furthermore,  the  data  acquisition  will  necessarily  cover  an  area 

instead  of  a  line  over  the  keel.  This  will  be  very  costly  and  may  not 
be  worthwhile. 

Errors  due  to  3)  need  to  be  further  investigated.  Treating  the 

ice  top  as  a  flat  plane  has  an  effect  of  pressing  the  ridge  down  into 
the  sea  water.  The  ice  thickness  probably  may  not  be  largely 


affected.  But  it  remains  to  to  be  seen  whether  a  topographic 
correction  is  necessary  before  a  2-D  inversion  is  attempted. 

Applying  the  constrained  smooth  inversion  to  the  scaled  sum  of 
the  in-phase  and  quadrature  components  of  the  real  data  appears  to 
be  an  attractive  simple  scheme.  If  the  phase  of  the  response  does  not 
ch  ge,  this  is  identical  to  scaling  the  in-phase  component  only.  The 
performance,  mainly  the  resolution,  of  this  scheme  will  deteriorate 
with  lowering  the  operation  frequency.  The  range  of  the  frequency, 
in  which  the  above  procedure  can  be  applied,  needs  to  be  defined  by 
further  investigation.  Effects  of  the  conductivities  of  the  sea  ice  and 
water  will  also  reduce  the  resolution  of  the  keel  geometry  expected 
from  the  assumption  of  transparent  ice  and  superconducting  sea 
water,  which  need  be  studied  in  the  future. 

All  the  computations  in  the  above  inversion  have  been  done  on 
a  IBM  RT  PC  computer,  which  has  a  comparable  speed  as  the 
Microvax.  An  typical  case  of  inversion  of  M  =  N  =  20  takes  about  4 
hours  of  the  CPU  time.  At  the  present  the  computer  code  is  not 
optimized  and  it  can  be  shown  that  its  optimization  may  sharply 
reduce  the  CPU  time  by  a  factor  of  5.  It  appears  that  in  flight  data 
interpretation  will  require  some  highly  specialized  computing 
equipment. 

In  all  our  examples  the  ice  keel  position  was  assumed  known. 
In  practice  the  position  of  the  ice  keel  may  be  obtained  from  the 
videocamera  record  of  the  associated  ridge.  If  this  is  not  possible  the 
AEM  data  itself  can  be  used  to  estimate  the  position  of  the  ice  keel 
prior  to  the  inversion.  To  do  this  the  influence  of  the  altitude  should 
first  be  removed  from  the  data  so  that  the  remaining  anomaly  will 
indicate  the  position  of  the  keel. 

4.4  Discussion 

The  2-D  inversion  scheme  can  successfully  recover  the  ice  keel 
information,  which  is  usually  lost  in  1-D  inversion.  Because  the 
electromagnetic  induction  is  a  diffusion  process,  sharp  edges  can  noi 


be  resolved  from  the  data  and  we  set  out  to  seek  a  smooth  ice  keel, 
which  retains  the  major  features  of  the  true  ice  keel.  Although  the 
scheme  is  designed  to  w'ork  for  data  collected  in  the  induction  limit 
with  transparent  ice  and  superconducting  sea  water,  it  can  be 
applied  to  low  frequency  data,  where  the  sum  of  the  in-phase  and 
quadrature  components  is  scaled  to  approximate  the  required  data. 
At  the  present  the  range  of  the  frequency,  in  which  this  procedure 
can  be  applied,  remains  unclear  and  needs  further  investigation. 
With  refinements  and  optimization  of  the  computer  program 
however,  there  is  no  doubt  that  this  algorithm  can  be  used  routinely 
to  process  field  data. 

5.0  The  Three-Dimensional  Keel 

Under  usual  circumstances  the  ice-water  interface  in  arctic 
seas  exhibits  the  same  irregular  topography  that  is  normally 

associated  with  land  forms.  In  areas  covered  by  young  ice  however, 

the  topographic  features  are  minimal  and  a  one-dimensional  data 
interpretation  technique  for  the  average  thickness  of  the  ice  sheet 
from  airborne  electromagnetic  data  (Kovacs,  et.  al.  ,  1987)  can  be 
used.  This  technique  fails  in  areas  of  even  moderate  ice-water 
interface  topography  (Kovacs,  et  al,  1987).  In  cases  where  the  axial 
length  of  the  keel  is  very  large  (compared  to  the  AEM  system 
altitude)  it  is  possible  to  remedy  this  situation  with  the  use  of 
interpretation  charts  that  allow  for  a  two-dimensional  description  of 
the  keel.  If  necessary,  a  2-D  data  inversion  outlined  above  may  be 
carried  out  to  delineate  ice  keels  in  more  details.  Because  the  tw'o- 
dimensional  data  interpretation  method  appears  to  produce 

adequately  accurate  results  it  appears  worthwhile  to  find  the 
minimal  keel  axis  length  that  is  necessary  for  its  proper  application. 
This  is  done  below  where  theoretical  data  for  three-dimensional  (3- 
D)  keels  is  presented.  In  all  14  different  keel  models  are  considered. 


5.1  The  keel  model 


Cartesian  coordinates  are  used  and  the  XY  plane  is  chosen  to 
coincide  with  the  flat  part  of  the  ice-water  interface  so  that  the  Z  axis 
points  vertically  down  into  the  sea  water.  The  ice  is  assumed  to  be 
transparent  to  electromagnetic  waves  and  sea  water  is  assumed  to  be 
a  perfect  electric  conductor.  Both  the  co-axial  (HX)  and  co-planar  (HZ) 
coil  systems  are  considered. 

The  3-D  Gaussian  keel  model  is  an  extension  of  the  2-D  model 
described  above.  Its  shape  is  given  by 


t(x,y)  —  AC  0.361  0.361  w“  (1) 

where 

t(x,y)  is  the  protruding  depth  of  the  ice  keel  into  the  sea  water, 
a  is  the  keel  drawdown  or  maximum  keel  thickness, 
w  is  the  keel  width  at  t  =  A/2  on  the  cross  section  y  =  0, 
s  is  the  keel  width  at  t  =  A/2  on  the  cross  section  x  =  0. 

Here  we  define  the  y  direction  as  the  strike  direction  of  the  ice 
keel  and  hence  designate  s  as  the  strike  length  of  the  keel.  An 
infinite  value  of  s  corresponds  to  a  2-D  ice  keel.  When  s  =  w  ,  the 
keel  is  a  body  of  revolution.  The  numerical  calculations  were  carried 
out  as  outlined  in  the  Appendix. 

5.2  Airborne  electromagnetic  profiles 

In  this  section,  we  examine  AEM  profiles  over  three  sets  of 
model  ice  keels.  Except  in  the  keel  area,  the  sea  ice  is  assumed  to 
have  a  uniform  thickness  of  5  meters.  The  keel  parameters  A,  w, 
and  s  are  varied  to  demonstrate  the  corresponding  changes  in  the 
AEM  system  response.  For  each  model  set,  the  keel  drawdown  A  and 
keel  width  w  are  fixed,  while  the  keel  strike  length  is  varied.  Model 
Set  1  represents  a  shallow  keel,  where  A  =  3m,  w  =  24m,  and  keel 
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strike  lengths  of  12,  24,  48,  96,  and  inf.  meters  respectively  (Figure 
20).  The  corresponding  theoretical  AEM  response  for  the  HX  and  HZ 
system  over  the  center  of  these  features  along  a  line  perpendicular 
to  the  keel  strike  is  shown  in  Figure  21,  where  curves  1,  2,  3,  4,  and 
5  correspond  to  s  =  12,  24,  48,  96,  and  inf.  meters  respectively.  Both 
the  transmitter  and  the  receiver,  which  are  separated  by  6.5m,  are 
"fl^wn"  25m  above  the  upper  icc  surface.  The  system  response  is 
plotted  at  a  point  located  mid-way  between  the  transmitter  and 
receiver,  and  is  expressed  in  parts  per  million  of  the  primary  field  at 
the  receiver.  Figure  21  (c)  shows  the  the  cross  section  of  the  ice 
directly  below  the  flight  line. 


•  It  is  evident  in  Figure  21(a)  that  as  the  strike  length  s  decreases 

the  HX  system  anomaly  becomes  smaller  but  the  general  shape  of  the 
anomaly  curve  changes  little.  If  these  system  anomalies  for  3-D  ice 
keels  are  interpreted  using  the  2-D  chart  given  in  our  previous 
0  report,  the  interpreted  cross  section  of  the  keel  will  be  smaller  than 

that  in  the  3-D  model.  The  significant  effect  here  is  a  reduction  of  the 
interpreted  keel  drawdown.  The  keel  width  however  is  less  affected 
since  the  anomaly  width  does  not  contract  significantly  with  the  keel 
strike. 


As  shown  in  Figure  21(a),  the  HX  system  response  for  a  keel 
strike  length  s  of  96m  (curve  4)  closely  resembles  the  response  for 
an  infinitely  long  keel  (curve  5).  In  fact  the  maximum  difference  in 
the  relative  amplitude  or  percent  anomaly  (c.f.  Becker  and  Liu)  is  of 
the  order  of  1.5%  and  would  result  in  an  under-estimation  of  the  keel 
depth  of  about  10%.  In  fact  the  relative  error  in  the  percent  anomaly 

size  is  also  of  this  order  and  suggests  that  a  10%  value  for  this 

quantity  be  used  as  a  threshhold  value  for  the  definition  of  a  two- 

dimensional  keel.  Thus,  If  the  maximum  relative  difference  in  the 
relative  anomaly  is  less  than  10%,  the  finite  strike  keel  can  be 

effectively  considered  as  two-dimensional,  and  hence  2-D 
interpretation  techniques  can  be  applied. 
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Using  this  criterion,  any  ice  keel  whose  strike  length  exceeds 
96m  can  be  effectively  considered  as  two-dimensional  for  the 
purpose  of  data  interpretation.  In  general,  the  ratio  of  the  strike 
length  to  the  system  height  over  the  ice-water  interface  must  be 
larger  than  96/30.  This  value  is  much  greater  than  the  value 
expected  from  the  footprint  for  the  HX  system  (see  Becker,  at  al., 
1987),  which  is  1.4.  At  the  first  sight,  this  may  seem  puzzling.  A 
careful  look  at  the  current  distribution  on  the  water  surface, 
however,  will  yield  the  explanation. 

The  current  pattern  for  a  horizontal-axis  transmitter  is  shown 
in  Figure  22  (from  Becker,  et  al.,  1987).  It  consists  of  two  circular 
current  rings  while  most  of  the  current  flows  in  a  direction  parallel 
to  the  strike  of  the  ice  keel.  For  a  keel  of  finite-strike,  this  major 
current  flow  is  distorted.  The  footprint  however  was  computed  on 
the  basis  of  the  unaltered  current  flow  so  that  it  is  not  surprising  to 

find  that  the  strike/height  ratio  needed  for  2-D  interpretation  is 

much  larger  than  the  footprint  of  the  HX  system. 

Now  let  us  look  at  the  HZ  system  response  and  see  how  it 
changes  as  the  keel  strike  length  is  varied.  Curves  1,  2,  3,  4,  and  5  in 
Figure  21(b)  correspond  to  s  =  12,  24,  48,  96,  and  inf.  meters 

respectively.  As  expected,  the  HZ  system  response  also  decreases 

when  the  strike  length  of  the  keel  is  reduced.  The  major  change  is 
the  decrease  of  the  anomaly  amplitude.  The  anomaly  width,  in 
contrast,  stays  almost  constant.  Nonetheless,  as  the  keel  strike  is 

further  decreased,  the  anomaly  shape  also  changes  and  the  lower 

part  of  the  anomaly  flattens  out. 

The  HZ  system  response  for  s  =  96m  is  also  close  to  that  for  the 
2-D  keel.  The  maximum  relative  difference  in  the  relative  anomaly  in 
this  case  is  about  12%  .  Thus  taking  10%  as  the  threshhold,  the  keel 
strike  should  be  somewhat  longer  than  96m  in  order  to  use  the  2-D 

interpretation.  Roughly  speaking,  we  may  again  take  the 

strike/height  ratio  to  be  approximately  96/30  =  3.2  for  the  2-D 
interpretation  to  be  valid.  This  is  the  same  as  that  for  the  HX  system 
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and  is  very  close  to  the  value  obtained  from  the  footprint  of  HZ 
system  which  is  3.7  .  As  a  matter  of  fact,  the  current  pattern  on  a  flat 
water  surface  for  a  vertical-axis  transmitter  consists  of  concentric 
circles  so  that  a  3-D  ice  keel  will  distort  this  pattern  as  much  as  a  2-D 
ice  keel  does. 

Now  let  us  increase  the  keel  drawdown  A  to  6m  and  keep  the 
keel  width  at  w  =  24m.  Figure  23  shows  ihe  deeper  ice  keels  with 
different  strike  lengths  (Model  Set  2).  The  corresponding  HX  and  HZ 
system  responses  are  shown  in  Figure  24  (a)  and  (b)  respectively. 
Now  both  the  HX  and  HZ  system  responses  are  much  larger  than 
those  for  Model  Set  1  due  to  the  larger  cross  section  of  these  keels. 
But  the  characteristics  of  these  anomalies  are  similar  to  those 
discussed  for  Model  Set  1.  Therefore,  the  analysis  for  Model  Set  1 
still  holds  in  this  case. 

Figure  25  shows  the  narrower  keels  that  make  up  Model  Set  3, 
where  A  -  3m  and  w  =  12m.  The  keel  strike,  again,  is  12,  24,  48,  96, 
and  inf.  meters  respectively.  Figure  26  (a)  and  (b)  displays  the  HX 
and  HZ  system  response  and  curves  1,  2,  3,  4,  and  5  correspond  to  s  = 
12,  24,  48,  96,  and  inf.  meters.  Figure  26  (c)  shows  the  cross  section 
of  the  keel  below  the  flight  line. 

Now  the  HX  system  response  is  much  smaller  than  those  shown 
in  Figure  21  (a)  because  of  the  smaller  size  of  the  ice  keels.  But  the 
relative  characteristics  remain  unaltered  and  the  analysis  for  Model 
Set  1  still  applies  here. 

The  HZ  system  anomaly  also  reduces  with  the  decrease  of  the 
strike  length  of  the  keel.  Furthermore,  it  is  observed  that  a  minor 
feature  begins  to  appear  at  the  center  of  the  anomaly  when  the  keel 
strike  is  small.  This  feature,  however,  may  not  be  real  as  it  may  be 
caused  by  numerical  errors  in  the  calculation.  At  present,  we  are 
unable  to  make  this  point  clear  due  to  computational  limitations. 


5.3  Discussion 


1)  The  strike  length  of  an  ice  keel  must  be  at  least  three  times  the 
flight  height  of  the  AEM  system  in  order  to  successfully  use  the  2-D 
techniques  to  interpret  both  HX  and  HZ  systems  response. 

2)  For  the  HX  system,  the  anomaly  amplitude  decreases  as  the  keel 
strike  shortens,  but  the  anomaly  width  changes  little. 

3)  For  the  HZ  system,  the  anomaly  amplitude  also  decreases  as  the 
keel  strike  contracts.  At  the  same  time,  the  lower  part  of  the 
anomaly  tends  to  become  flat. 
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Fig.l  HX  and  HZ  systems  over  sea  ice 
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Fig. 9  Interpretation  Chart  for  HZ  system 
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Fig.  10  Field  data.  Line  F6L3. 
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Figure  12  (a)  Synthetic  data  for  the  HX  system  and  inversion  results.  Model 
(b)  Synthetic  data  for  the  HZ  system  and  inversion  results.  Model  1. 
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Figure  17  (a)  Height  of  the  system  over  the  ice/water  interface  during  the 

flight,  (b)  Synthetic  data  and  inversion  results  for  the  MX  system.  Model  3.  ( c ) 
Synthetic  data  and  inversion  results  for  the  HZ  system.  Model  3. 
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Figure  18  Synthetic  data  at  2500Hz  and  inversion  results  for  the  HX 
system.  Model  4.  The  data  shown  are  the  scaled  sum  of  the  in-phase 
and  quadrature  components  of  the  secondary  magnetic  field. 


t—  cn  J 
X  ^  . 
CD  _ 


60  9 

distance:  (m 


Q_  CD 

Q-  g- 


^  ,  o 
m  o 
x  o 


t— 

Q-  CD. 

LJ 

CD 

CD  ■ 


WATER 


__  scaled  data 
response  of  inverted  keel 


60  90 

DISTANCE  (M) 


drill-hole 

inverted 


Figure  19  Scaled  in-phase  dita  at  16,290Hz  and  the  inversion  results 
for  the  HZ  system.  The  data  were  collected  at  the  Prudhoe  Bay  by 
Geotech  Ltd.  in  1985  on  line  F6L3. 


Figure  21(a)  HX  system  response. 
Curves  1,  2,  3,  4,  and  5  corres¬ 
pond  to  s  =  12,  24,  48,  96,  and  ini 
meters  respectively. 
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Figure  21(b)  HZ  system  response. 
Curves  I,  2,  3,  4,  and  5  corres¬ 
pond  to  s  =  12,  24,  48,  96,  and  in 
meters  respectively. 
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Figure  21(c)  Cross  section  of  the  ice 
below  the  flight  line. 


Y(M) 


50 


o 


ri 


o 


CM 


O 

CM 


O 
Q  " 


O 
CM  - 
I 


O 


OJ 


* 

/ 

/ 

i 

i 

k 

V 

> 

> 


'  '  N  V  \ 

-  «•  *>  N  V  \  \ 

►  »  '  v  N  \  \ 

-  ^  V  S  \  \  \ 

x  \  \ 
«■  -  **  X  \  \ 
*•  -  «*  v  \  \  \ 

✓  »  ^  s  \  \  \ 

/  -  *  v  '  \  \ 

*  '  **  '  '  \  \ 

*  i  1  »  «  *  i 

V  V  -  '  t  i  ] 

>  '  -  '  '  i  I 

'>-•*■'//  / 

'  -  ■*  ^  /  /  / 


o 

CO 

ro 

I 


-36.0 


- T - j 

-24.0  -12.0 


'  1  ♦  » 
*  i  i  t 
\  l  *  / 
\  \  \  t 
\  \  / 
\  \  / 

ill 


/It 

//it 

//It 
/  /  M 
/  f  »  * 
/  t  *  » 
/ft* 


0.0 

XIM) 


/  /  //>>--. 

/  /  /  y  /  -  -  .  . 

/  /  ///■•--» 

/ 

( 

///✓---'> 

/  ///»---.> 

I  j  /  /  ^  >  v 

/  /  f  '  -  V  v  v 

|  i  *  »  ’  »  »  * 

i  ^  \  4  i  i 

\  \  \  *  ~  *  *  < 

\  \ 

\  \  \  \  v  -  r 

\  \  V  V  v  •>  '  -  » 

S  \  \  \  s  K  » 


12.0 


24.0 


Figure  22  Surficial  currents  for  a  horizontal-axis  transmitter. 

which  is  30  meters  above  the  center  (0,  0). 
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Figure  24(a)  HX  system  response 
Curves  1,  2,  3,  4  correspond  to 
s  =  24,  48,  96,  and  inf.  meters 
respectively. 
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Figure  24(b)  HZ  system  respons^ 
Curves  1,  2,  3,  4  correspond  to 
s  =  24,  48,  96,  and  inf.  meters 
respectively. 
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Figure  24(c)  Cross  section  of  the  ic< 
below  the  flight  line.  * 
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Appendix  A 


Computational  Methods 


A.l  The  General  Case 

Consider  an  alternating  magnetic  dipole  (  current  loop  )  source 
T.  located  in  free  space  as  shown  in  Figure  Al.  Its  orientation  is 
arbitrary  and  it  is  positioned  over  a  homogeneous,  perfectly 
conductive  medium  with  three-dimensional  surface  relief.  Due  to  this 
source,  electric  currents  are  induced  on  the  surface  of  the  medium 
and  they  give  rise  to  the  secondary  magnetic  field  Hs  in  free  space.  It 
is  our  objective  to  calculate  this  quantity  at  any  point  above  the 
surface  S. 

Here  the  height  of  the  source  above  the  medium  is  assumed  to 
be  small  compared  to  the  wavelength.  The  observation  point  is  also 
assumed  to  be  close  to  the  source  so  that  the  electromagnetic  field  is 
quasi-static.  In  free  space,  V  x  H$  =  0  ,  and  we  may  relate  H  s  to  a 
scalar  magnetic  potential  (J) 

Hs=-Vx0  (1) 

We  also  have  V  •  H$  =  0  and  correspondingly 

V24>  =  0  (2) 


Since  the  lower  medium  is  assumed  to  be  infinitely  conductive,  the 
normal  component  of  the  total  magnetic  field  must  vanish  on  the 
surface  S,  i.e. 

Hpn  +  Hsn  =  0  on  S  (3) 

Here,  Hpn  and  Hsn  are  the  normal  components  of  the  secondary 
and  primary  magnetic  fields  respectively.  Hence, 
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_ 


-is=-Hjs  =  h 


ipn's 


The  Laplace  equation  (2)  and  the  boundary  condition  (4) 
constitute  the  Neumann  boundary  value  problem.  That  is,  given  the 
normal  derivative  of  the  potential  on  a  surface  S,  we  wish  to 
calculate  the  potential  itself  in  the  free  space.  Once  <{)  is  found,  Hs 
may  be  calculated  from  equation  (1). 


The  solution  of  the  Neumann  problem  outside  a  closed  surface 
can  be  expressed  as  the  potential  of  a  surface  charge  layer  (Graham,  # 

1980) 


<>(o)  = 


Tpo 


(5) 


where  ^(P)  is  a  fictitious  charge  density  function,  rop  is  the  distance 
between  points  P  and  O  (  see  Figure  1  ).  The  charge  density  satisfies 
a  Fredholm  integral  equation  of  the  second  kind 


S(M)  =  - 


1  3d)  | 
2k  dn  s 


cos(rPM,N)  ± 

TFsi 


on  S 


(6) 


where  (  rPM,  N  )  is  the  angle  between  rPM  (  the  vector  connecting  P  to 
M  )  and  N  (  the  unit  normal  vector  at  M  ).  In  our  problem,  the 
surface  S  extends  to  infinity  and  equations  (5)  and  (6)  are  still  valid. 


To  solve  the  integral  equation  (6),  we  use  the  successive 
approximation  method  (  Mikhlin,  1964  ).  The  initial  solution  (first 
iteration)  is  assumed  to  be  the  first  term  on  the  right  hand  side  of 
equation  (6),  i.e. 


£(M)  =  - 


2k  9n  s 


2k 


We  then  use  this  value  in  the  integral  in  the  equation  to  compute  an 
improved  value  of  %(M)  and  so  on.  Once  the  charge  density  is  known. 


the  secondary  magnetic  field  can  be  calculated  from  the  following 
equation,  which  is  obtained  by  combining  equations  (1)  and  (5) 

Hs=  |  5Mr,ods 

-h  r'°  (  7  ) 

Here  rpo  is  the  vector  connecting  P  to  O. 

In  the  above  context,  the  electrodynamic  problem  is  reduced  to 
a  potential  problem  under  the  quasi-static  field  assumption.  The 
charge  distribution  on  the  surface  of  the  conductive  medium  is 
computed  first.  The  secondary  magnetic  field  in  the  free  space  is 
then  found  by  summing  up  the  contributions  from  the  individual 
charges.  This  is  analogous  to  the  integral  equation  approach  for 
solving  electromagnetic  scattering  problem  (Parry  and  Ward,  1970), 
where  the  equivalent  electric  and  magnetic  currents  are  first  sought; 
the  electromagnetic  field  are  then  obtained  by  integrating  the 
contributions  from  the  current  distribution. 


In  the  special  case  where  the  surface  S  forms  a  plane,  the 
solution  given  by  equation  (7)  for  an  oscillating  magnetic  dipole 
source  is  identical  to  that  obtained  by  the  method  of  images  (Jackson, 
1975).  To  illustrate  our  computational  method,  we  now  show  this  to 
be  true  for  a  vertical  magnetic  dipole  source.  The  coordinate  system 
is  chosen  such  that  the  z  axis  is  pointing  vertically  downward.  The 
XOY  plane  is  on  the  surface  of  the  conductor  which  occupies  the  half 
space  z  >  0  .  The  dipole  source  is  h  meters  above  the  plane  on  the  z 
axis  and  points  in  the  positive  z  direction. 


Since  rpM  is  perpendicular  to  N  (cf.  Figure  1),  equation  (6) 
reduces  to 


4(M)  =  -  J-  |*  |  =  -  J_  H 
2k  9n  *  271 


pn 


(8) 


On  the  conductor  surface,  the  normal  component  of  the  primary 
magnetic  field  is 
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H  =  1  2h2-x2-y2 

*  47T  (X2+y2+h-)5 


Substituting  equations  (8)  and  (9)  into  (5),  we  obtain 


<>(x,  y,  z)  = 


- S1*.'.-  /’  0) - dx'dy’ 

[(x-x')-+(y-y'):+z2]1/2 


(z<0) 


0(0,0,  z)=  -L 

8it- 


2h~-xl2-y'~ _ dx'dy' 

(x’2+y'2+h2)5/2  (x’2+y'2+z2)1/2 


Let  x'  =  Rcos9  ,  y'  =  Rsin9  ,  dx'dy'  =  Rd9dR,  then 


0(0,  0,  z)= 

8  TL~ 


0  JO 


(R2+h2)572  (R2+z2)172 


Furthermore,  let  r~  =  R"-  +  h-  ,  RdR  =  rdr  ,  then 


0(0,0,z)=J-  - ihyrL^dr 

4KJh  ^(^-h'+z2)172 


The  above  integration  can  be  carried  out  and  the  result  is  as  follows 


0(0,  0,  z)=  -J - 1— 

4k  (h-z)2 


(z<0) 


The  vertical  component  of  the  secondary  magnetic  field  on  z  axis  is 


14,(0, 0,  z)  =  -^-=J-— !— 
oz  2k  (h-z)3 


(10) 


and  as  required  is  identical  to  the  image  field. 
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A. 2  The  2-D  keel 

In  the  2-D  case,  the  geometry  of  the  ice-water  interface  does 
not  change  in  the  strike  direction  (  y-direction  here  ).  The  relief  of 
the  interface  is  only  a  function  of  x  ,  i.e.  h(x.y)  =  h(x).  In  this  case, 
the  potential  of  the  scattered  magnetic  field  is  given  by, 

^(x'.y') Vi  +[dh(x')/dx'D  ,  ,,  , 

■  ■  ■  —  —  ax  ay 

V(x-x')2+(y-y')2+(z-h(x'))2  (H) 


0(x,  y,  z) 


-If! 


and  the  surface  charge  density  £(x,  y)  satisfies 

-lx  _J_H  f-:  1  Pf-:'  v"i  lx-x,)dh(x’)/dx'  ~  (h(x)-h(x')) 

[(x-x')2+(y-y')2+(h(x)-h(x'))‘]2 

.  rilwWilWirydx.dy. 

l+(dh(x)/dx)‘  >p 

Here  Hpn(x.y)  is  the  normal  component  of  the  primary  magnetic  field 
at  the  ice-water  interface. 

Notice  that  in  equations  11  and  12  the  integral  with  regard 
to  y'  is  a  convolution.  Taking  the  Fourier  transform  of  both  sides  we 
get 


0<X  k  z)  =  J__  0(x,  y.  z)  e‘  lk>'y  dy 


:  2.f_^  ^(x’,ky)  Vl+[dh(x')/dx']2  K0Cplkyt)  dx’ 


(13) 


and 


1  1 

q(x,  k  )  = - H  tx,ky) - J  £(x',k  )  f(x,x',ky)  d> 


»  *pn' . y>  _ 

2k  2  k 


where  the  kernel  of  the  integration  is 


(K) 
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c,  ■  r  \  -  -,rl+(dh(x')/dx')2-|i-  (x-x')dh(x')/dx'  -  (h(x)-h(x’)) 

I(X,  X  ,  Ky)  -  4.  - Jl - - - 

l+(dh(x)/dx)2  f(x-x')2  +  (h(x)-h(x'))2]2 

Jkyl  Ki  (p'  Iky  I )  (^15) 

In  the  above  equations 

\r 

y  =  angular  wave  number  in  the  y  direction 

p = V  (x-xV+ (z  -  h(x')  y 
P'=V(x-x’)2+(h(x)-h(x’))2 

Ko(  )  =  modified  Bessel  function  of  the  zeroth  order.  Second  kind. 

K-i(  )  =  modified  Bessel  function  of  the  first  order.  Second  kind. 

We  have  thus  simplified  the  problem  by  decomposing  a  two 
dimensional  integral  equation  into  a  number  of  integral  equations  in 
one  dimension. 


For  the  case  where  the  source  is  a  horizontal  or  vertical 
magnetic  dipole,  the  normal  component  of  the  primary  magnetic  field 
can  be  analytically  transformed  into  the  wave-number  domain  (i.e. 
Hpn(x,  ky)  )  as  follows. 


The  normal  component  of  the  primary  magnetic  field  at  the 
ice-water  interface  is 


Hpn(x.y)  =  Hp(x,  y)-n(x,  y) 


(16) 


where  n(x,  y)  is  the  outward  normal  at  point  (x,  y)  on  the  interface. 
It  is  given  by 


n 


h’fx) 


,  0. 


Vl  +[h'(x)]2  ’  ’  Vl 


[h'(x)r 


where  h’(x)  =  dh(x)/dx.  Thus 


(17) 
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Hpn(x,  y)  = 


h:(x) 


Vi 


H 


1 


+  [  h’  ( x)  ] 


P* 


V 


H 


1  +  [h'(x)]‘ 


P' 


(18) 


(a)  Horizontal  magnetic  dipole  source  located  at  (xs.  ys=0,  zs) 


Hpx(x,y)  = 


2(x  -  -  y2  -  (h(x)  -  zy 


Hpz(x,  y)  = 


4?tr5 


3(x  -  xe)  (h(x)  -  zj 


4nr5 


Here 


r  =  V  (x-Xs)“  +  y2  +  (h(x)  -  z,p' 


Substituting  the  above  equations  into  (18)  and  taking  the  Fourier 
transform  give 


H  (x,kv)  = - f  [  h'(x)(x-Xs)2  -  (x-Xs)(h(x)-zs)]  Ik yl2  K0  (p!kyl) 

2np‘Vl+  fh'(x)J"  '• 

-  ([(hfxj-Zj)2  -  (x-x j) 2  ]  h'(x)  +  ICx-XjaHW-zJ)  K,  (plkyl)  | 


(19) 


w  h  e  r  e 


p  =  V  (x-Xj)2  +  (h(x)  -  zj2 


As  ky  — >  0,  the  asymptotic  forms  of  the  modified  Bessel  functions 
are 

K0(p|kyl)  -  In  plk  1 


K  i  (p|k.  I) 


plkyf 


Therefore 


-  1 


HnJx,k=0)~ -  r - 

Pn  >  4  f  1 

2?tpVl+  fh(x>] " 


=•  {[(h(x)-zv)2  -  (x-x/  J  h'(x)  +  2(x-Xs)(h(x)-zs) } 


(20) 
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(b)  Vertical  magnetic  dipole  source  located  at  (xs,  ys=0,  zs) 


y)  ~ 


3(x  -  xj  (h (x)  -  zs) 


4%r5 


H  (x,k  )  = - --F  ==•  (  [-(h(x)-Zj)2  +  h'(x)  (x-x^hCxJ-z,)]  Ik  A2  K0  (p!k  I) 

2^pVl+[h'(x)]2  ^ 


[(h(x)-z,)  -  (x 


-Xj)2  -  2h'(x)  (x-XjXhCxJ-Zs)]  —  Kt  (plkyl)  | 


(21) 


H  (x,k  =0)  = - 7- 7^======  {(h(x)-Zj)2  -  (x-Xj)2  -  2h'(x)  (x-x^hM-zJ ) 


27rp4Vl+[h'(x)]2 


(22) 


Note  here  that  the  kernel  f(x,  x',  ky)  is  independent  of  the 
source  and  the  receiver  positions.  Therefore  it  may  be  calculated 
once  for  each  ky  value  and  stored  for  the  computation  of  a  complete 
profile  of  AEM  system  response.  This  is  quite  economical  but 
impossible  in  the  3-D  case  where  the  matrix  is  too  large  to  store.  At  x 
=  x',  the  kernel  has  a  singularity.  But  this  presents  no  difficulty  for 
the  numerical  computation  because  it  is  integrable  in  the  sense  of 
Cauchy  principal  value. 


The  integral  equation  (14)  can  be  solved  for  q(x,  ky)  using 
successive  approximation  method  identical  to  the  one  suggested 
above  for  the  general  case.  This  needs  to  be  done  at  a  number  of 
positive  ky  harmonics  (including  ky  =0).  The  values  of  c,(x,  ky)  at  the 


negative 


ky  harmonics  may  be  easily  obtained  by  its  property 


of 


symmetry.  As  usual,  the  sampling  in  the  ky  space  is  done  on  a 


logarithmic  scale. 
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Once  the  charge  density  is  known,  the  x-  and  z-  components  of 
the  scattered  magnetic  field  may  be  directly  computed  from  the 
following  equations 


Hsx(x,  ky,  z)  =  - 


30(x,  kv,  z) 


=  2/^  (  x  -x'  )  ^(x’,k>)  VT+idh(x')/dx'] ‘  ^-K^p'ky!)  dx’ 


and 


Now  that  the  scattered  magnetic  field  is  obtained  in  the  w-ave- 
number  domain  its  mversf*  Fourier  transform  will  yield  the  desired 
result  in  the  space  domain.  But  prior  to  performing  the  inverse 
Fourier  transform,  me  field  values  at  the  logarithmically-spaced 
points  in  the  kv  space  need  to  be  interpolated  for  uniformly  spaced 
values.  This  is  accomplished  using  cubic  spline  interpolation. 


The  above  algorithm  has  been  successfully  implemented  and 
its  speed  is  more  than  20  times  faster  than  that  of  the  general  3-D 
algorithm.  The  computation  of  a  20  point  profile  of  the  AEM  system 
response  takes  about  only  10  seconds  CPU  time  on  the  IBM  3090. 


A. 3  Computational  Check 


A  scale  model  experiment  was  also  performed  to  check  the 
numerical  solution.  Aluminum  was  used  to  simulate  the  infinitely 
conductive  medium  at  a  scale  of  1:250.  A  model  airborne 
electromagnetic  system  was  built  at  the  same  scale  and  was  "  flown 
at  a  field  height  of  10m.  The  system  consisted  of  a  coplanar, 
vertical  axis  transmitter  and  receiver,  which  w'ere  operated  at  6  kHz. 
Details  of  the  model  are  shown  in  Figure  A2  which  also  exhibits  the 
traversed  feature. 
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The  cross  section  of  the  indented  surface  is  a  Gaussian  curve 
that  simulates  a  smooth  ice  keel.  Its  relief  is  given  by 

t(x)  =  Aexp(-  a!) 

2xz  (25) 

and 


x  =  0.425W 

Here 

x  =  distance  from  the  keel  center  line 
t  =  keel  thickness 

A  =  drawdown  or  maximum  keel  thickness 
W  =  keel  width  at  half  drawdown 


The  shape  of  the  keel  does  not  change  along  its  strike  direction.  For 
this  scale  model,  A  and  W  were  taken  to  be  3.4m  and  21m 
respectively.  We  chose  this  type  of  surface  because  its  simulation  on 
the  computer  is  simple  as  there  are  only  two  input  parameters. 
Furthermore,  it  is  easy  to  adjust  these  two  parameters  to  simulate 
any  real  sea  ice  keel.  The  measurements  and  numerical  calculation 
results  (two  iterations)  for  this  model  are  displayed  in  Figure  A3  , 
which  shows  an  excellent  agreement  between  the  numerical  and 
experimental  data. 


0 


filR 


Fig.Al  The  Neumann  Problem 


Distance  (m) 


Fig.  A3  Comparison  of  Measurements  and  Computation 
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•  -point  Gauss'  quadrature 
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data  w/0 . 330239 ,0.312347,0. 260611 ,0 . 180648,0.081274/ 
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The  i^pilization  of  downhole  current  sources  in  resistiv  ity  mapping  increases  the  resolution  for 
detecting  and  delineating  subsurface  features  The  effects  of  near  surface  inhoniogeneities  are 
immensely  reduced  as  shown  by  Asch  and  Morrison  (1988)  Being  sensitise  to  changes  in  resistivi¬ 
ties,  the  surveys  with  downhole  sources  are  well  suited  for  monitoring  surface  processes  such  as 
injection  or  leakage  of  contaminants  from  a  waste  site,  steam  flooding  for  enhanced  oil  recovery,  or 
production  of  geothermal  reservoirs 

In  most  of  these  applications,  the  holes  are  steel-cased  and  the  saving  distorts  the  current  flow 

■jl-  y  ■ 

in  the  medium.  Holladay  and  West  ( 1984 1  ha8r  shown  the  surla  e  resistivity  surveys  are  strongly 

r 

'-"aftected  bv  casines.  Also,  Kaualukauj  et  al  ildtsoi  showed  that  it  the  casine  itself  is  used  as  an 
electrode,  the  results  are  unpredictable  because  the  current  leaves  the  pipe  irregularly  due  to  the  vari¬ 
ability  of  the  contact  resistance  between  the  pipe  and  the  formation 

To  study  the  casing  effects  in  more  detail,  we  have  recently  formulated  the  problem  for  a  point 
source  of  current,  either  inside  or  outside  the  pipe,  on  the  avis  of  a  finite  length  metal  pipe  in  a  con¬ 
ductive  half  space.  The  first  part  of  this  study  (Schenkel  and  Morrison.  1987)  showed  that  only  the 
region  very  near  the  pipe  exerts  any  substantial  influence  on  the  potential  fields  for  a  point  source 
100  casing  diameters  beyond  the  end  of  the  pipe  (Figure  I  )  For  3  591  or  less  field  distortion,  the 
surface  measurements  must  not  be  closer  than  1/2  the  pipe's  length.  In  cross-hole  surveys,  the 
affected  area  is  greatly  reduced;  near  the  pipe’s  end,  measurements  as  close  as  1/6  the  length  of  the 
pipe  for  a  591  or  less  distortion.  If  the  pipe-source  separation  is  sufficiently  large,  then  the  resistiv  ity 
survey  can  be  corrected  for  the  casing  effects  (assuming  that  the  target  are  not  too  close  to  the  pi  pc ) . 
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Thi>  siudy  showed  that  cross-hole  and  hole-surface  studies  may  be  conducted  with  very  little 
effects;  even  when  the  pips'  exerts  some  influence,  e  g  .  when  the  current  source  is  close  to  the  end. 
time  monitoring  experiments  could  be  carried  out  with  very  little  reduction  in  the  signal  strength 
Further,  this  studs  lesealed  the  intriguing  possibility  that  segments  of  pipe,  separated  b>  insulated 
couplings,  could  be  used  as  electrodes. 

Proposed  Work 

The  above  studs  was  only  developed  for  a  vers  simple  case,  i  e..  a  pipe  in  a  homogeneous 
half-space  A  more  complex  model  is  required  to  simulate  field  situations.  Several  aims  are  pro¬ 
posed  to  create  a  realistic  model  and  to  evaluate  field  measurements  for  downhole  sources  in  steel 
cased  wells.  These  objectives  are: 

li  To  determine  the  effects  of  contact  resistance  between  the  pipe  and  the  host  medium.  Contact 
resistance  is  used  to  describe  the  resistance  of  the  pipe-medium  contact  If  there  is  a  large 
contact  resistance,  the  results  will  completely  change  The  currents  in  the  pipe  will  only  leak 
out  of  the  pipe  in  areas  where  the  pipe  has  made  a  good  contact  with  the  host  medium,  thus 
completely  changing  the  potential  fields  around  the  pipe.  The  contact  resistance  may  be  found 
by  assuming  a  very  thin  layer  between  the  host  and  pipe  for  each  segment  (Figure  2).  An 
equivalent  layer  can  be  used  to  represent  the  effects  of  the  contact  resistance  for  each  sege- 
nient  The  calculated  potential  field  is  obtained  from  the  integral  equation  solution  of  the  pipe 
variables  With  the  additional  layer  included  to  the  model,  the  effects  of  pipe  coating,  corro¬ 
sion.  and  cement  on  the  potential  fields  can  be  evaluated 

2)  To  use  insulated  pipe  segements  as  downhole  sources  and  receivers.  Downhole  current  elec¬ 
trodes  can  be  created  by  energizing  isolated  segements.  Likewise,  insulated  segements  could 
be  used  as  potential  electrodes.  By  insulating  several  segments  in  the  well  (Figure  3),  multiple 
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downhole  source  locations  could  he  used  to  image  a  target  in  a  hole-surlace  survey  A  pole- 
pole  contiguration  can  be  acheived  by  attaching  to  different  segments  current  and  potential 
electrodes  The  isolated  segments  in  the  well  can  be  used  for  an  AC  vertical  electric  dipole 
By  attaching  to  adjacent  segments  a  positive  and  negative  AC  source,  a  grounded  electric 
dipole  can  be  produced  to  studs  F.M  properties  of  the  medium  If  an  additional  well  is  drilled 
with  multiple  isolated  segments,  then  cross-hole  DC  and  AC  measurements  can  be  acquired. 
Thus,  cross-hole  DC  tomography  (Daily  and  Morels,  1988  and  Shiina  and  Satlo,  1988)  could 
be  used  to  reconstruct  an  image  of  a  target  between  the  two  wells. 

3i  To  determine  the  interaction  between  the  source,  pipe,  and  the  anomalous  bod)  (Figure  4i  To 
monitor  the  changes  tn  the  bods',  the  effects  of  the  metal  pipe-body  interaction  must  be  investi¬ 
gated  The  extent  of  (fits  interaction  will  greatly  depend  on  the  source  separation  from  the 
pipe,  the  distance  between  the  pipe  and  the  body,  and  the  conductivity  of  the  body.  The  deter¬ 
mination  of  the  limiting  values  of  these  parameters  in  which  the  body  and  pipe  has  very  little 
coupling  will  be  the  main  objective.  For  this  situation,  the  body  can  be  modeled  alone  saving 
computational  time. 

4|  To  invert  for  the  geometric  parameters  for  a  plume-like  body  An  integral  equation  solution  of 
an  ellipsoidal  model  will  be  used  to  represent  the  plume.  The  parameters  of  the  three  axial 
lengths  will  be  obtained  by  a  non-linear  least  squares  inversion  which  will  make  use  of  the 
integral  equation  solution  Sensitivity  analysis  and  minimum  spatial  coverage  will  be  evaluated 
for  various  acquisition  array  configurations. 

Computer  models  will  be  required  to  investigate  the  above  proposed  tasks.  The  current  algo¬ 
rithm  is  flexible  so  that  variable  source  locations,  segment  lengths,  and  segment  conductiviiies  can  be 
used  The  development  of  an  algorithm  which  includes  an  anomalous  body  and  an  outer  layer  is 
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needed  The  current  integral  equation  solution  can  he  extented  to  include  an  additional  laser  and  a 
circular  disk  The  circular  disk  is  a  first  order  approximation  to  a  plume  and  will  give  an  estimate  of 
the  pipe-bods  interaction 

Various  lengths  and  separations  of  the  insulated  segments  will  be  investigated  to  determine 
when  a  point  approximation  of  the  segments  can  be  used.  The  outer  iayer  will  be  used  to  evaluate 
the  role  of  contact  resistance  on  the  field  distortion  and  mas  be  included  «>  the  pipe-bod>  coupling 
phenomenon.  The  spatial  separation  of  the  source,  pipe,  and  bod>  to  decouple  the  pipe  and  the  bodv 
w  ill  be  studied. 

Field  test  and  model  verification  are  also  required.  The  test  would  be  conducted  at  U.C 
Berkelev.  Richmond  Field  Station  where  there  exist  several  plastic-cased  holes.  Four  additional  100 
foot  holes  are  also  needed  Two  of  the  welts  will  be  composed  of  alternating  steel  and  fiberglass 
segments  Another  would  be  cased  with  steel  with  the  last  20  feet  perforated.  The  last  well  needs  to 
be  plastic  cased  and  perforated  at  the  bottom  A  surface  grid  and/or  radial  arrays  of  potential  elec¬ 
trodes  would  be  installed  ever  the  area  Cross-hole  and  hole-surface  measurements  would  be  com¬ 
pared  to  calculated  fields  of  various  pipe  models.  An  injection  of  the  salt  water  in  both  the  steel- 
cased  and  plastic-cased  wells  would  be  measured  for  hole-surface,  cross-hole,  and  pole-pole 
configurations.  These  data  would  be  forward  modeled  and  inverted  to  determine  the  geometry  of  the 
plume  Lastly,  other  configurations  for  which  models  have  been  published  can  be  field  tested  with 
this  field  setup. 
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